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Abstract 

We perform a detailed and illustrative study of the production of keV sterile neu¬ 
trino Dark Matter (DM) by decays of singlet scalars in the early Universe. In the 
current study we focus on providing a clear and general overview of this produc¬ 
tion mechanism. For the first time we study all regimes possible on the level of 
momentum distribution functions, which we obtain by solving a system of Boltz¬ 
mann equations. These quantities contain the full information about the production 
process, which allows us to not only track the evolution of the DM generation but 
to also take into account all bounds related to the spectrum, such as constraints 
from structure formation or from avoiding too much dark radiation. In particu¬ 
lar we show that this simple production mechanism can, depending on the regime, 
lead to strongly non-thermal DM spectra which may even feature more than one 
peak in the momentum distribution. These cases could have particularly interesting 
consequences for cosmological structure formation, as their analysis requires more 
refined tools than the simplistic estimate using the free-streaming horizon. Here 
we present the mechanism including all concepts and subtleties involved, for now 
using the assumption that the effective number of relativistic degrees of freedom is 
constant during DM production, which is applicable in a significant fraction of the 
parameter space. This allows us to derive analytical results to back up our detailed 
numerical computations, thus leading to the most comprehensive picture of keV 
sterile neutrino DM production by singlet scalar decays that exists up to now. 
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1 Introduction 


Despite great advances in our understanding, our Universe holds mysteries that are yet to 
be resolved. One of the biggest questions is about the identity of the so-called Dark Matter 
(DM) which - in terms of cosmic energy density - outweighs ordinary matter by a factor 
of about five [lj[2|. While historically (motivated by, e.g., supersymmetry) our best guess 
for DM was a Weakly Interacting Massive Particle (WIMP) with a mass of a few hundred 
GeV and roughly weak interaction strength, by now we have unfortunately not seen a clear 
signal of such a particle. Even worse, attempts for direct detection disfavour big parts of 
the parameter space that was deemed “natural” |3[[6]. While WIMPs cannot be excluded, 
we are nevertheless at a point where we should seriously think of alternatives |7|. 


After all, there are several possibilities left for what DM could be. Alternative ideas range 
from very light scalars such as axions [ 8 | over non-standard fermions in supersymmetry 


(such as gravitinos 9| or axinos |!10|) up to very heavy exotic options like WIMPzillas 11 


In this work, we concentrate on a candidate motivated by neutrino physics, a sterile 
neutrino u s of (typically) a mass of a few keV. While our natural guess would be for the 
sterile neutrino mass to be much larger, it is in reality not bound and there exist several 
consistent theoretical frameworks in which sterile neutrinos can be very light [12,13 


Sterile neutrino DM has originally been proposed by Dodelson and Widrow (DW) [14] 
who suggested that, although sterile neutrinos would not have interactions strong enough 
to keep them in thermal equilibrium in the early Universe, they could nevertheless be 
produced gradually by the thermal plasma due to their admixtures to active neutrinos. 
Although the DW mechanism was at that time consistent with all bounds 15 , by now 
we know that it is excluded by structure formation: it produces too hot DM 16 . Taking 


the sterile neutrino mass to larger values is not possible due to the X-ray bound on the 
DM decay v s —>• v + 7 , where v is any active neutrino (see Ref. [l7| for a comprehensive 
discussion and Ref. 18] for probably the most recent collection of limits)^] 


Several other production mechanisms for keV sterile neutrino DM have been proposed. If a 
primordial lepton asymmetry is present in the early Universe, the active-sterile transitions 
could be resonantly enhanced, as proposed by Shi and Fuller (SF) [28]. This mechanism 


Tn 2014 two groups have independently reported a tentative signal of an X-ray line at 3.5 keV 


19 


20 


which would, if interpreted as originating from sterile neutrino decay, indicate a particle with a mass of 
7.1 keV and an active-sterile mixing angle 9 of roughly sin 2 (29) ~ 7 ■ 10 -11 . However, up to this date 
there is still a discussion on-going in the literature about whether or not this line is in fact a real signal, 
see Refs. 21-27 . Given that the signal, if it exists, is at the moment still not at a statistically highly 


significant level, we do not take any side here but instead suggest to wait until more data is collected, to 
hopefully either strongly confirm or clearly refute a line signal. 
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produces a relatively cold DM component in addition to the thermal DW part, which 


leads to a colder spectrum in accordance with all bounds 17,29,30 - however, note that 


also this mechanism is in some tension with data if the 3.5 keV line is taken seriously 31 


In theories with an extended gauge group, the sterile neutrinos could equilibrate and be 
produced via generic freeze-out, if the resulting overabundance is corrected by a sufficient 
amount of entropy production 32,33 - however this scenario, even though not fully 


excluded, is on quite general grounds threatened by big bang nucleosynthesis 34 


An alternative is non-thermal production of keV sterile neutrinos via particle decays. 
This possibility is discussed extensively in the literature, with the decaying particle in 


most cases being a scalar: variants range from inflaton decay 35,36 over a general scalar 
singlet that could itself either freeze-out 37,38 or freeze-in 39-41 to the case where 


the scalar is electrically charged 42 


the production from pion decays 43 


More general possibilities exist as well, for example 
from Dirac fermions 


44 , or from light vector 


bosons 45,46. A benefit of this mechanism is that the velocity spectrum of the keV 


neutrinos produced by decays is quite generally known to have a tendency to be colder 
than that from other mechanisms |3l][38,39i|47j48]. 

While several cases of decay production are discussed, the treatments available in the 
literature involve some crude estimates and approximations. Even though a non-thermal 
DM spectrum could have interesting and unexpected consequences for cosmological struc¬ 
ture formation, many references work on the level of rate equations and only estimate the 
spectrum, if at all. Instead, given that the distribution function is the most decisive 
quantity and that it can be computed from Erst principles with reasonable effort, its de¬ 
termination should be put on more solid grounds. This is partially done for the production 
of keV neutrinos by a general singlet scalar entering thermal equilibrium 38 , however, 


only in the approximation where the effective number of relativistic degrees of freedom 
g* is constant. While this approximation is not usable in all of the parameter space, it is 
nevertheless a good approximation for large enough production temperatures and, after 
all, without this assumption it would not be possible to obtain analytical estimates. Yet, 
ultimately a fully numerical treatment will be necessary. 

In this paper, we will start this endeavour by giving a comprehensive discussion of the 
production of keV sterile neutrino DM via the decays of general singlet scalars, which 
themselves are produced via a Higgs portal. As shown in Refs. 38-40 , depending on 


the coupling there exist different regimes where the scalar itself could e.g. be a WIMP 
or a feebly interacting massive particle and decay in or out of thermal equilibrium. We 
will generalise the previous treatments and derive the full set of approximate formulas 
for all regimes possible, but again under the assumption of g * = const. We furthermore 
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perform a fully numerical study of the allowed regions in the parameter space and obtain 
distribution functions for all relevant parameter points, which in principle allows us to 
determine all relevant DM properties for a given point in the parameter space. We 
determine the regions where the correct DM abundance is obtained for a given sterile 
neutrino mass, thereby taking into account all relevant bounds. 

Structure formation properties of DM candidates can be estimated using the so-called free- 
streaming horizon Aps- With this quantity, we show that different estimates are possible 
which may lead to quite different results. We clearly illustrate that Aps, in spite of being 
a standard estimator, can lead to inconsistent results depending on the approximations 
used to calculate it. 


The current paper serves as an illustration of the general principles and its purpose is to 
give a clear overview of the relatively complicated and subtle details involved. As such, 
some aspects lie beyond the scope of the current work. For example, we neglect active- 
sterile mixing and thus completely disregard DW production, which in reality cannot be 
switched off for non-zero active-sterile mixing (and which on the contrary is desired for the 
related X-ray signal), in favour of analytical results. Our approximation is still good for 
very small active-sterile mixing, but it is nevertheless too simplified and will be dropped 
(as well as the assumption g * =const.) in a future purely numerical investigation of all 
possible cases 49 . Also the full set of implications for cosmological structure formation 


cannot be obtained using Apg, so that a numerical computation of the structure formation 
properties is needed. While we could in principle add this to the current work, it would 
lead away from our main point and furthermore prolong the paper even more, so that 
we have decided to decouple this study, too 50 . Our goal is to ultimately deliver a fully 


comprehensive study of scalar decay production of keV sterile neutrino DM, in order to set 
the stage to discriminate it from alternative mechanisms by future data. Current bounds 
indicate that this may be possible in the not-too-far future 31 , so that this endeavour 


should be pursued in particular if the 3.5 keV line signal survives. 


This paper is organised as follows. We start with an illustrative discussion in Sec. [2j 
which is supposed to give an overview of our considerations and results at a relatively 
non-technical level. We then introduce the main equations to be solved in Sec. i in 
Sec. [4j we present analytical approximations for all cases where they can be obtained. 
After a discussion of the aspects relevant for cosmic structure formation, Sec. [5j we finally 
present the numerical results of our study in Sec. [6j along with a discussion of all bounds 
and of the validity of our considerations. We conclude in Sec. [7} Throughout the paper 
we try to keep the discussion on a minimally technical level, however, all technical details 
relevant for an inclined reader are exposed in Appendices A and B. 
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2 Illustrative discussion of the setting and overview 
of the results 


This section mainly serves the purpose of describing in simple terms what we have done in 
our study. Before entering the technical details, we give an overview not only of the setting 
we are working in but also of some illustrative results. This will hopefully make the paper 
more accessible and prevent the reader from getting lost in the technical details. 

Our basic setting consists of a real singlet scalar field S which must somehow couple to 
the right-handed neutrino fields N. The most generic coupling doing this job is a Yukawa 
term | SN C N with coupling strength y which, if the scalar develops a non-zero vacuum 
expectation value (VEV) (S), leads to a Majorana mass mjv = y{S), where we have 
assumed only one generation of right-handed neutrinos for simplicity. Thus, the minimal 
set of ingredients we need in addition to the standard model (SM) consists of exactly 
these fields. Our minimal Lagrangian is 


£ — £sm + 


iN$N + ^ (d^S) (d^S) - | SN C N + h.c. 


Scalar Y C-v 


( 1 ) 


which is very similar to what had been used previously 38-40 . Here, C v denotes the 


part of the Lagragian giving mass to active neutrinos. The details of this mass generation 
are in fact not very important for our DM production mechanism, as is the number of 
fermion generations. However, in the current work, we make the additional assumption of 
vanishing active-sterile mixing. In most settings, there will be couplings between active 
and sterile neutrinos that in reality cannot be switched off. But, since in this paper we 
want to present analytical results wherever possible, we take this mixing to be exactly zero 
and in turn have no DM production from the DW mechanism 0 The scalar potential Vji Ca i ar 
takes the most general form compatible with an assumed global Z 4 -symmetryj^] 


Scalar = ~ H ~ ^S 2 3 + X H H)* + ^ S ' 4 + 2A (rf H) S 2 


( 2 ) 


2 Given that the astrophysical X-ray bounds push the active-sterile mixing down to tiny values 17 18 


this approximation is not necessarily bad. However, an additional contribution from the DW production 
can modify some of the results obtained here, which is why we will drop this assumption in future works 


and turn to a purely numerical treatment 49 50 . 

3 For possible issues related to the breaking oFthis symmetry by a non-vanishing VEV (5), see 


39 


and 


references therein. Note also that giving up the Z 4 symmetry allows for extra terms in the Lagrangian, 
resulting in more processes that can contribute to equilibrating the scalar. In this paper, we stick to the 
symmetry in order to obtain a minimal parameter space for our exploratory analysis. Considerations of 


taking into account terms linear and cubic in S in the potential can be found in 38 . 
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This potential could easily lead to a VEV (S) ~ GeV — TeV, thus motivating a relatively 
small Yukawa coupling y ~ 10 -9 -10~ 4 5 in order for the mass of the sterile neutrino to be 
in the keV-range. 


Not only can the above Yukawa coupling lead to a sterile neutrino mass, but it will also be 
responsible for sterile neutrino DM production. Several processes could be thought of: the 
leading contributions are the reactions SS O hh and S —y NN, the first of which couples 
the scalar field S to the thermal plasma and the second of which produces sterile neutrinos 
from the decay of S. In principle, also processes like SS —» NN would be possible, but the 
corresponding rate is proportional to y A which is negligible for a sufficiently small Yukawa 
coupling y. We also neglect the inverse reaction NN —> S which is suppressed due to the 
heavily suppressed phase space originating in the kinematics of any 2-to-l process. 


Let us add that, in general, there will be mixing between the scalar S and the SM Higgs 
field after electroweak symmetry breaking, which we have completely ignored. However, 
given that this mixing is proportional to the generally small Higgs portal coupling A and 
that it is also suppressed if the singlet scalar is considerably heavier than the Higgs, which 
is the case that we are investigating here (cf. Sec. |4]) , taking into account the mixing would 
not at all change our results. Note that this would change if one used singlet scalar masses 
very close to or below the Higgs mass |40|, which is why this simplifying assumption must 
be dropped if we are to extend our considerations to lighter singlets. On the other hand, in 
that limit we would need to drop further assumptions used in this work (in particular the 
assumption that g* is constant), so that it makes sense to postpone these considerations 
to future work [49]. 


With these assumptions, it is clear that every scalar present in the early Universe will 
either decay into two sterile neutrinos or undergo pairwise annihilation into pairs of Higgs 
bosons. Depending on the exact values of the couplings, there are different regimes 
possible. For example, if the Higgs portal coupling A is small enough and the initial 
abundanc^] of the scalars is zero, then the scalar will never enter thermal equilibrium 
but it will only be produced occasionally from the plasma. In a more modern language, 
this mechanism would be called freeze-in 5l] 52], and the corresponding particle would 
be called a feebly interacting massive particle (FIMP). In this regime, the annihilation 
into Higgs bosons can be neglected since its reaction rate will be suppressed by the square 
of the (tiny) scalar density. Hence, the frozen-in abundance of scalars will ultimately be 
translated into a relic abundance of sterile neutrinos with a particle number just twice 


4 We usually use the term abundance to denote a particle number density. Depending on the context, 

the term relic abundance will be used for the particle number density or the corresponding energy density 
after the process discussed is complete. 
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the one of the scalars at freeze-in (in a co-moving volume) - irrespective of the size of 
the Yukawa coupling y as long as it is large enough that all scalars have decayed by 
now. Another example would be the case where the scalar couples to the Higgs strongly 
enough to be in thermal equilibrium. Then the argument of doubling the number of 
particles would still be valid if the decay width of the scalar is sufficiently small for the 
decay to become effective only after freeze-out. However, if the decay proceeds already 
while the scalar is in equilibrium, further contributions will be present making the whole 
picture more complicated. We will discuss all cases in detail in Sec. [4j where we present 
analytical estimates for the momentum distribution functions / and yields Y wherever 
possible. 

We later on turn to a numerical computation of the DM production. To arrive at the final 
DM relic abundance, we need to solve a system of Boltzmann equations to compute the 
momentum distribution function f(p,t) of the DM particle. Ultimately this distribution 
function contains all relevant information: one can e.g. use it to compute the final DM 
relic abundance, but it also encodes the information about the evolution of the momentum 
spectrum with time, i.e., how many particles exist per momentum interval at any given 
temperature of the Universe. This allows to not only track the course of DM generation 
in detail, but it furthermore gives information about the velocities of the DM particles at 
the epochs relevant for the formation of structures (i.e., galaxies and galaxy clusters) in 
space. 

To give a snapshot of the results to be presented, we show in Fig. [l] a plot of the lines 
of correct DM abundance for several different sterile neutrino masses, in a plane of the 
parameters Chp and Cr which, as we will explain later, are nothing else than an effective 
Higgs portal coupling and an effective decay width in convenient units. We have aug¬ 
mented the plot by some example evolutions of the DM production and the underlying 
distribution functions for some specific parameter points marked in the central figure, in 
order to illustrate what is behind our computations. The purpose of Fig. [T| is to give 
a graphical illustration of what will be presented in this paper. All the plots displayed 
will be discussed in great detail in Sec. [6j where also the terminology, colour-codes, and 
labelling will be carefully outlined so that, while advancing with the paper, the reader 
will ultimately be enabled to get a full understanding of Fig. [T[ 


6 


Hr)=nir)!T 3 



Figure 1: Lines of correct abundance with example distributions and evolutions. 



























































3 The kinetic equations in the early Universe 


The fundamental embodiment of every particle species in the early Universe is their distri¬ 
bution function f(p, t) in momentum space. This quantity does contain all the information 
we need to deduce the cosmological impact of the species under consideration, so that the 
determination of / is paramount. Due to isotropy, we will always assume that the distri¬ 
bution functions / only depend on the moduli of the associated 3-momenta. To compute 
a distribution function, we need to solve the corresponding Boltzmann equation: 


i[/] = c[/]. 


(3) 


The left-hand side of this equation contains the so-called Liouville operator L: 


f _ d_ _ d_ 
dt Hp dp ’ 


(4) 


where p is the modulus of the particle’s 3-momentum and H is the Hubble function. The 
collision term C [/] on the right-hand side can be interpreted as a source term, encoding 
the interaction of the species of interest with itself and with the other particle species 
present in the plasma. This term contains all the information about the processes which 
contribute to the production of the species under consideration, and which accordingly 
shape the resulting momentum distribution function. Collision terms can be relatively 


lengthy, which is why we report the explicit form of C [/] only in Appendix A In this 
section, however, we prefer to give a more intuitive explanation of how to obtain the 
distribution functions of the sterile neutrinos. 


Since DM production happens in the very early Universe, we only need to consider the 
era of radiation dominance. During this epoch, the Hubble function can be written as 
H (T) = T 2 /M 0 , where T is the temperature of the plasma and M 0 is a function that 
implicitly depends on time via the evolution of the degrees of freedom gp. 

M o = 7 = 7 - 3 %* _1/2 x 10 18 GeV. (5) 

Introducing dimensionless variables, 


x=p/T and r = ms/T , 


(6) 




the Liouville operator from Eq. Q can be recast into the following form: 

f dr dT_5_ H \ Q_ 

“dT d tdr X {l-^£V9~*J dx ' 


( 7 ) 


Throughout this work, we will stick to the assumption that the numbers of effective degrees 
of freedom ( g*,g*s ) are constant until the production of sterile neutrinos is completed , as 
done e.g. in [38,[39]. This assumption is absolutely necessary to obtain analytical results 
and, in fact, it is not a too bad approximation in a large fraction of the relevant parameter 
space, as we will illustrate in Sec. |6.5| Nevertheless we will drop it in later works where 


we are going to present more realistic and purely numerical studies 49,50 


If the number of effective degrees of freedom does not change during the period of interest, 
the advantage of the variables x and r becomes obvious. Accordingly, the dynamics of 
the scalar and the sterile neutrino are given by the following set of equations: 


dfs 

dr 


dT d^ 
dr dT 


(c> 


ih- 


■ss 


+ C S 


SS^-hh 


+ c s 


S^-NN) ) 


( 8 ) 


df N _ dT d t N 
dr dr dT C ' S ^ NN ' 


(9) 


In Eqs. (j8j) and (|9|, the upper indices on the collision terms mark the species the kinetic 
equation of which is governed by this term, while the subscripts describe the actual 
process. Thus, Cg S ^. hh describes the depletion of scalars due to annihilations into pairs 
of Higgses while C% h ^ ss describes the reverse process. In turn, C<j^ NN describes the 
depletion of scalars due to decays and Cg_^. NN encodes the creation of sterile neutrinos 
from the decays of scalars. Note that the collision terms contain information about the 
kinematics, too, so that Cg^ NN and C§^ NN differ by more than just a sign. For a detailed 
derivation and explicit expressions, see Appendix [A] 

Since we approximate g * as constant during the time of production and since the matter- 
radiation equality only takes place at temperatures of O (leV), i.e., long after DM pro¬ 
duction is completed, we consistently use the time-temperature relation ^ = —HT in 
Eq. ([9]). Using this as well as the explicit form of the collision terms as derived in Ap- 
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pendix |A.1[ we find the kinetic equation for the scalar: 


dfs (x, r ) 
dr 


\J x 2 +' 


-^-C H p exp (y—\/x 2 + r 2 j J 7 (or, r) — C r r 2 f s (x, r ) 


( 10 ) 


1 


— — Chp/s (x, r) 27 t J dxx 2 J dcos9fs(x,r)Q(x,cos9,r) 
o -i 

= Q (x, r) - V (x, r)f s -n (x, r) X r [f s ] fs , 


where we have defined 


C H p exp (-Vx 2 + r 2 ) T (x, r) 
Q [x, r) = - 


47Ta/x 2 + r 2 


V (x,r ) = 
7 Z (x, r ) = 


C r r 2 


\Jx 2 + r 2 ’ 
^hp 

47r\/a; 2 + r 2 ’ 


I r [/g] = 2tt dx x 2 / dcos 9fs (x, r) Q (x, cos 9, r). 


-l 


( 11 ) 

( 12 ) 

(13) 

(14) 


For the explicit forms of the kinematic functions T and Q and the definition of a max , see 
Appendix [A} 


In Eqs. © to dl3|, we have used two important dimensionless auxiliary quantities: 


the effective decay width: 

the effective (squared) Higgs portal. 


C 

C 


= Mo, 

1 ms ms ’ 
_ Mp A 2 
ms 167T 3 


(15) 


In the remainder of this work it will turn out convenient to use these quantities to span 
the parameter space of our setting. Note that, during the DM production process, we 
assume M 0 to be constant by virtue of Eq. ([5]). Hence the interpretations of Cpp as an 
effective Higgs portal coupling and of C p as an effective decay width are appropriate, in 
the sense that the dependence of Mo on g * does not change their behaviour and hence 
they indeed play practically the same roles as the fundamental Lagrangian parameters 
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behind them. For simplicity, we may often refer to Chp as Higgs portal and to Cp as decay 
width. 


Turning to the sterile neutrino distribution function, which is our main quantity of interest, 
we can use the assumption of neglecting the back-reaction NN —y S to find a very simple 


for details, one can derive a very intuitive master equation for the distribution function 
of the sterile neutrino in terms of that of the scalar: 


form for the corresponding kinetic equation. Using Eqs. (pi) and rt9J), see Appendix A.2 


f n (x, r) 



V 

dr' 2Cp — 5 - 
x z 



x 


V x 2 + r' 2 


fs (£, r '), 


(16) 


with f m in = ||x — r , 2 /(4a;)||. It is this equation that will be absolutely central for us: 
once we have managed to determine the distribution function fs of the scalar from first 
principles, we only need to plug it into Eq. (16) to determine the distribution function 
f N of the sterile neutrino. As to be expected, the distribution function fs of the scalar 
will look differently depending on which regime we are looking at (e.g., the scalar being 
an early decaying WIMP compared to a late one). This will be translated by Eq. (16) 
into different resulting sterile neutrino distribution functions. Note that, transforming 
the momentum variable into an energy, Eq. (16) perfectly coincides with 35, Eq. ( 8 )] 
and 38, Eq. (10)]. Note also that the master equation for the sterile neutrino distribution 
is decoupled from the one for the scalar only because of our assumption that the reaction 
NN —>■ S is negligible. 


4 Analytical results for the distribution functions and 
relic abundance 

In this section, we present the limiting cases that can be treated analytically. Some of 
the results can be found in the literature, see in particular Ref. [38], while others are 
completely new. Basically there exist two main cases, the scalar itself can either enter 
thermal equilibrium in the early Universe and thus act similarly to a WIMP or it can be 
only very feebly interacting, thus undergoing freeze-in like a FIMP: 

1. The FIMP-regime-. This limiting case is characterised by a small Higgs portal Chp 
and an (almost) arbitrary decay constant Cp. The scalars that are produced via 
freeze-in subsequently decay into sterile neutrinos. However, for the abundance 
itself the time of the decay is not very relevant as long as it happens before matter- 
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radiation equality, which is why the exact value of Cp does not matter much in this 
regime]^] It will, however, play a role for the spectrum itself, cf. Sec. [6} 

2. The WIMP-regime: For large enough Higgs portals Chp the scalar will thermalise, 
i.e., it enters thermal equilibrium and any information about its initial abundance 
is lost since the abundance will always be the thermal one. The particle remains 
in equilibrium until the interaction rate for the process hh •<->• SS drops below the 
expansion rate of the Universe and the scalar decouples from the thermal bath. 
During the course of these events, depending on the exact value of the decay width 
Cp, the scalar can decay into sterile neutrinos at various stages: 


(a) In-equilibrium decay. If the decay width Cp is large enough and if the Higgs 
portal Chp is large, too, the sterile neutrinos are produced already very early 
from the decays of equilibrated scalars. The sterile neutrinos arising from the 
decays of the few residual scalars that are present after freeze-out can then 
practically be neglected since the large Higgs portal guarantees a small relic 
abundance of frozen-out scalars. 


(b) Out-of-equilibrium decay. The scalar couples to SM particles strongly enough 
to enter thermal equilibrium. However, if the decay width is sufficiently small, 
the production of sterile neutrinos during the time of equilibrium is negligible 
and it is sufficient to only take into account the late decays of the scalars after 
their freeze-out. In this regime, the scalar itself acts as a DM-like species, but 
it is unstable and decays before it can significantly contribute to the energy 
budget of the Universe. 


(c) Intermediate regime: For intermediate values of the decay width Cp, neither 


5 Note that, for this regime, we will always assume a zero initial abundance for the scalars (and trivially 
for the sterile neutrinos). If there was a non-negligible initial abundance, this would change our results 
since the primordial scalars would then add to the ones produced via freeze-in. However, given that 
there is no reason for such an initial abundance to be present and that we do not see much value in 
speculating how it could possibly have been produced, we stick to the conservative viewpoint and only 
produce scalars from the freeze-in mechanism itself. 

On the other hand, one could argue that sterile neutrinos and/or singlet scalar fields could quite generically 
couple to the inflaton field (see Refs. 35 36 for examples concerning the former case). In such scenarios, 
assuming that inflation is the correct theory in the first place, our scenario might be modified considerably. 
However, such couplings are only compulsory if the SM gauge group and all other low-energy symmetries 
are the only ones up to very high scales and even then, from a model building point of view, there exist 
various possibilities to strongly suppress certain couplings, e.g. by locating the various fields on different 
branes. While our considerations could in principle be extended to include this point, this would add 
further complications of which it is however unclear whether they exist, or not. We thus stick to the most 
minimal setting and disregard any primordial abundance of sterile neutrinos and/or the singlet scalar, as 
well as the related coupling to the inflaton field. 
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of the above limiting cases is applicable and it is a priori not clear how well 
the combination of both of them can correctly describe the situation. On 
the other hand, this case can be of particular interest since it results into 
a distribution function with two intrinsic momentum scales. This may open 
up interesting possibilities to tackle the well-known small scale problems of 
cosmological structure formation 53-56 . We will therefore treat an explicit 


example for this case numerically in Sec. 5.2 


We will now discuss these different regimes one by one. In all cases, we will take the limit 
£ = rrih/ms —>• 0, cf. Appendix A.l Even for £ « 0.3, the net effect on the distribution 
function is of the order of a few per rnille. We want to keep the mass ms of the scalar far 
from the electroweak (EW) phase transition to ensure that scalars S are only produced by 
Higgs bosons. For masses ms < i>ew, a new range of interaction becomes important like 
the production of EW gauge bosons from scalars S, SS yy VV. This scenario is discussed 
in 40 on the level of abundances, and we will take these modifications into account on a 


further more technical study where the focus will be put on even more realistic numerical 


results 49 . 


4.1 The FIMP-regime 


For sufficiently small Higgs portals Crp (he., for A <C 10 6 38,39]), the scalar never enters 
equilibrium and hence one can - for vanishing initial abundance - neglect in Eq. (10) the 


term 1Z (r, x) T r [/$] which scales as /J. The remaining kinetic equation is an ordinary 
differential equation which can be solved analytically. For a vanishing initial abundance 
of the scalar, the resulting distribution function is given by: 


fs {r, x) = C 


HP 


dp pK\ (p) 


exp 




X‘ 




X z 


9 V^+^ ( P + V /^T 


X A 


e rVr*+x* 1 r + ^ r 2 + x 2 


C r /2 


(17) 

where K\ is the first modified Bessel function of second kind, cf. Eq. (A-10). Plugging 


Eq. ( |T7| ) into Eq. (16) in order to get an analytical result for the distribution function of 
the sterile neutrino is not very instructive, since there is no simple form of that expression. 
However, the abundance of the sterile neutrino can be computed analytically for late times, 
r —* oo. Setting Cp to zero corresponds to a stable scalar. With this choice, Eq. © 
can be integrated rather easily and one obtains the (hypothetical) relic abundance of the 
stable scalar. Since all frozen-in scalars will however ultimately decay into two sterile 
neutrinos for a non-vanishing decay width C r , the abundance of sterile neutrinos will then 
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just be twice the abundance of the would-be-stable scalar 39 . The result for the yield 


Y = n/s, with n and s the particle number and entropy densities, respectively, is given 
by 

Chp 


Yn M oo) = 


135 


64vr 2 g* (T prod ) ' 

Here, T prod denotes the temperature at the time of production^ 


(18) 


4.2 The WIMP-regime 


In-equilibrium-decay If both the Higgs portal Chp and the decay width Cr are large, 
sterile neutrinos are efficiently produced from the decays of scalars already while being 
in equilibrium. This case is covered at length in 38 . However, our analytical expression 


differs by a constant factor, which is why we will sketch the most important steps needed 
to derive the result. If the scalar is in equilibrium, we know its distribution function 
exactly. Accordingly, the authors of 138j use a Bose-Einstein (BE) distribution to capture 
the quantum nature of the scalar. Since the whole set of equations governing the dynamics 
was however derived using the Maxwellian approximation in the Boltzmann equation, one 
might wonder whether it would be more consistent to again use a Maxwell-Boltzmann 
(MB) distribution for the scalar. We will exemplify both cases in order to illustrate that 
the difference between them is irrelevant. 


Plugging the BE and MB distributions into (A-20) yields: 


f N ( x , r) = 


8C r f dzxy/z- Hog ( 1 _ e 1 _ 3;z ) (BE) 


8C r f dz x-s/z — le xz = 


-xz _ e (x/A^A)) 


2^/x 


c xz ' yfz r - 1 (MB) 


(19) 

where we have introduced the variable z r = r 2 /( Ax 2 ) + 1 for convenience Q Integrating 
/at (x, r) over d 3 x and again taking the limit r —>• oo allows to calculate the yield for late 


6 Of course the time of production is subject to some ambiguities in its definition. Both freeze- 
in/freeze-out of the scalar and its subsequent decay are continuous processes, the time scales of which 
are determined by Chp and Cr- It is hence convenient to define the production time as the point when 
the abundance of sterile neutrinos has passed some threshold fraction of the final abundance, which we 
take to be 95%. 

'Note that, as we had already mentioned, we have neglected the mixing between the two physical 
scalars, which is a very good approximation in our case. However, in order to simplify the comparison of 
our results to the ones obtained in Ref. 


to the results from that reference. 


38 , it is of course necessary to apply the same approximation 
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times: 


Y n (r —> oo) = 


^< 5 >^ < BE > 

&— (MB) 


135 


( 20 ) 


9 * ("Jprod) 

Both results only differ by a factor of ( (5) ~ 1.0369, which justifies the use of either 
distribution. Our result in the BE case is however larger by a factor of 5/2 compared to 
the one reported in 38 . While one may easily forget powers of two in these computations, 


we could not trace any step where a factor of 5 could possibly be introduced, making us 
confident that our results are correct. 


Out-of-equilibrium decay This limiting case describes a scenario where the scalar is 
in equilibrium and ultimately freezes out. The decay width C F is so small that practically 
no sterile neutrinos are produced before the scalar decouples from the plasma. Only 
after the scalar has frozen-out, the production of the sterile neutrinos sets in. As in 38 


we make the approximation that the scalar has a thermal distribution until it freezes out 


instantaneously at r = r FO , In this case the kinetic equation, Eq. (10), can be solved: 


fs (x, r > r FO ) = / eq (x, r FO ) 


r + \/r 2 + x 2 
Wo + V r FO + x2 


C T x 2 /2 


3 -C r (rVx 2 +r 2 -r FO y/x 2 +r 2 0 J/2 


( 21 ) 

where f eq (x, r F o) is the equilibrium distribution of S at freeze-out. It could again be 
taken to be BE or, more consistently, MB. In the case of BE, our expression coincides 


with 38, Eq. (43)]. The final abundance of sterile neutrinos can in this limiting case again 
be calculated from doubling the abundance of scalars at freeze-out. Given as a function 
of rpo, the expression for the yield is 


Y n {r —t oo) 


45 

47r 4 g* (T prod ) 


OO 



dee 


V e2 - r lo 

e e — 5 


rpo 


/_ 45rp 0 /i 2 (r FO ) 
V 47T 4 g* (T prod ) 


for MB 


( 22 ) 


with <5 = 1 (5 = 0) for the BE (MB) case. Also here the numerical difference between 
both versions is fairly small for realistic values of r F o. 


Intermediate regime If C F is in an intermediate regime, neither the production be¬ 
fore nor the one after freeze-out completely dominates the sterile neutrino distribution 
function, such that no simple analytical treatment is possible. In this instant we have to 
rely on a purely numerical treatment. We will present a sample case for this regime in 
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Sec. [6j Nevertheless this intermediate regime may be of special interest as it features two 
intrinsic scales. 


4.3 Calculating the relic abundance 


So far we have shown formulae to compute the yield at r —y oo. To convert this into 
the commonly quoted closure parameters, we first have to multiply the yield by today’s 
entropy density of so = 2891.2 cm -3 57 and then compare the DM energy density today 


to the critical density of the Universe. We can write 


Draw h 2 — 


?7ijvhjv(r —>■ oo)so 

Pcrit/ h 2 


(23) 


where p cr a/h 2 = 1.054 x 10 2 MeV cm 3 |57|. With this conversion formula at hand, it is 


straightforward to transform the analytical estimates for the yield Y/v hi Secs. |4.1| and |4.2 
into expressions for the DM relic abundance, which can then be compared to the observed 
la range obtained by the Planck collaboration, D DM h 2 = 0.1188 ± 0.0010 |2|. 


5 Aspects of structure formation 

Up to now, we have only been concerned with the DM relic abundance, i.e., the amount 
of DM in the Universe. However, for a viable DM candidate it is not only important to be 
sufficiently abundant in the Universe, but it must also allow for structures such as galaxies 
to form. Clearly, this imposes constraints on the momentum distribution function f(p,t) 
of the DM candidate under consideration: it must not be hot, i.e., it is not allowed to be 
too relativistic at the time of structure formation when galaxies form (more precisely, the 
allowed amount of hot DM is at most a tiny 1% of all DM in the Universe 

The form of the distribution function /(p, t ), i.e. the spectrum of the DM candidate, can be 
constrained from the observed matter distribution in the cosmos: the evolution of spatial 
inhomogeneities depends on the spectrum of the DM particles and therefore it ultimately 
constrains the DM production mechanism. Since it is computationally impossible to run a 
simulation of structure formation for any possible distribution function, a commonly used 
indicator that can be calculated easily for a given f(p,t ) is the so-called free-streaming 
horizon Afs- This quantity describes the average distance a DM particle would have 
travelled after production without collisions and not subject to gravitational clustering. 
In fact this quantity can serve as a good estimator for a length below which the formation 
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of structures is heavily suppressed. For that reason, the free-streaming horizon Afs is 
commonly used in the literature to discard DM models with a spectrum that is too hot 
to explain the filamentary structure of the large scale matter distribution in the Universe 
by preventing galaxy-sized objects of being formed. 


5.1 Treatment of the free-streaming horizon 


In this section, we will show that the free-streaming horizon itself is subject to substantial 
uncertainties in its definition, which will make clear why we later present some of our 
results in two different versions. Moreover we will demonstrate that even our simple one- 
component DM model can produce highly non-thermal momentum distribution functions, 
which may even feature two intrinsic moment urn/velocity scales. In such a case, the 
average momentum (and hence the free-streaming horizon) cannot be expected to lead to 
sensible conclusions. 

The free-streaming horizon is defined by [l6|: 


-Vs 


[ T ‘ (v (T)) At 
a(T) dT 


(24) 


where T pro( j is the temperature at which production can be seen as complete (in the sense 
discussed before) and T 0 is today’s temperature. Here, (v (T)) is the average velocity of 
the sterile neutrinos that can be calculated from the distribution function: 


(v(T)) = 


f dx , x3 f N (x,r ) 

0 \J x 2 +r 2 (m N /m s ) 2 

f 0 °° dx x 2 f N ( X , r) 


(25) 


From Eq. (25) it becomes clear that (v) converges to unity, i.e., to the speed of light, as 


long as fjy is concentrated around values of x — p/T y 7 r 2 m 2 N / m 2 s = m N /T , just as 
expected. Since fa (x) does not change after the production process is complete, the factor 
of r 2 m 2 N /m 2 s in the square root increases. Note also that, once r 2 m 2 N /m 2 s is greater than 


the value(s) of x around which / is concentrated, Eq. (25) converges to the non-rclativistic 
expression, (v) —> T/m n (x) = ( p) / m at. 

For our numerics, we follow the approximations usually found in the analytical ap¬ 
proach 16 , namely we will assume the Universe to be completely radiation dominated 
until some temperature T eq (“matter-radiation equality”), where the Universe switches to 
being completely matter dominated. The last epoch of vacuum dominance is irrelevant: 
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even though this period dominates the past of the Universe on an absolute time scale (the 
matter-vacuum equality being located at a redshift of z — ~ 2.2, corresponding 

to a time of roughly 3 x 10 9 a 1601), the velocities of the sterile neutrinos in this epoch are 
so tiny that the resulting contribution to the free-streaming horizon is negligible. Com¬ 
monly, the evolution of the degrees of freedom (d.o.f.) is implemented by an additional 
dilution factor of ,U/ 3 = [(?* 5 (T pro d)/( 7 *s(To)] 1 / 3 ~ (106.75/3.36) 1//3 ~ 3.17, by which Afs 
is rescaled to account for entropy dilution (cf. (38,39]). Note that, although we approxi¬ 
mate ~ const, during DM production, we nevertheless have to take into account the 
entropy dilution until today, since there is no justification for the above assumption to be 
valid through the entire history of the Universe. Departing from the above approximation 
also during DM production would modify the dilution factor £ 1//3 , but not too drastically 
because of the presence of the third root. We will later on perform an estimate of the 
validity of our approximation, cf. Sec. 6.5 Note that, however, in this formalism the 


dependence on g*s (T pro d) is quite mild, such that it is safe to use the SM number of d.o.f. 
g*s = 106.75 even though the new particles contribute as well to some degree, depending 
on how strongly suppressed their true distribution functions are compared to a thermal 
one. 

There is, however, an issue with the analytical approach. The above treatment does not 
capture the physical fact that, in reality, the evolution of the d.o.f. also enters in the time- 


temperature relation inside the integral in Eq. (24), but this can be taken into account 


rather easily in a numerical evaluation of the integral. We therefore compute Aps in a 
second (numerical) version, in order to take the full evolution of the d.o.f. into account. 
Thereby, our numerics uses a set of fitted analytical formulae for the evolution of the 
d.o.f. [6l]. For more technical details on this numerical integration, see Appendix [B] 

In order to compare to the results already present in the literature, we will follow both 
approaches in parallel. To get an idea about whether the DM can be classified as cold, 
warm, or hot for a certain set of parameters, we take Afs = 0.1 Mpc to mark the boundary 
between hot and warm DM. This choice is relatively common in the literature (see, e.g., 
Refs. 


62-651), and it corresponds to the size of a typical dwarf satellite galaxy, which 


yields a sensible physical motivation. The boundary between warm and cold DM in turn 
is smooth but it is clear that Aps should be considerably smaller in this case, so that 
we choose Afs = 0.01 Mpc, i.e., one order of magnitude smaller than for the hot/warm 
boundary. Of course there is some arbitrariness involved in these choices, but given 
that the free-streaming horizon in itself can only yield an indication, the actual error 
introduced is not as serious as it may seem at first sight. In general the free-streaming 
horizon can only serve as an order-of-magnitude estimate, and it clearly should not be used 
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to prematurely discard unclear results. Ultimately, only more advanced computations of 
structure formation can assess whether scenarios with borderline free-streaming horizons 
should be discarded, or not 50 . For recent developments beyond the commonly used 


free-streaming approach, see e.g. 66 


5.2 Free-streaming horizon failing — an explicit example featur¬ 
ing twin peaks 


As mentioned before, even our simple one-component DM model can feature a distribution 
function with two intrinsic scales, namely in the case where the relic abundance of sterile 
neutrinos is produced from the decay of equilibrated scalars and the decay of frozen- 
out scalars in comparable amounts. We present here the exemplary case of Chp = 10 4 
and Cr = 1CD 3 ' 5 , which yields the correct relic abundance for a (relatively large) sterile 
neutrino mass of about 73keV. Fixing the scalar mass to be 1 TeV and the number of 
degrees of freedom at high tempertures to g* = 106.5, this would correspond to values of 
the Lagrangian couplings of (y, A) = (4.7 x 1CT 9 , 8.3 x 1CT 5 ). 

Fig. [2] shows the distribution function of the sterile neutrino for different values of the 
time parameter r. One can clearly see how early production from the plasma populates 
the lower comoving momenta (dubbed Thermal part - although the resulting distribution 
may not be perfectly thermal), while the late contributions mainly originate from the 
decay of the frozen-out scalars (Decay part). The inset in the plot shows the enlarged 
region between the two peaks. 


It is obvious that in this case the average momentum does not at all capture the char¬ 
acteristics of the distribution function. According to our estimates using Afs and the 
average value (x) = (p) /T ~ 16.6, this point in the parameter space corresponds to a sce¬ 
nario where the sterile neutrinos are on the borderline between being hot and being warm 
(cf. Sec. [6]). However, in fact the low momentum (“thermal”) part with (x)\ ow ~ 2.5 con¬ 
tributes practically as strongly as the high momentum (“decay”) part with (x)hi g h ~ 35.7, 
where in both cases we have approximated the respective peaks with the individual results 
from Eqs. (19) and (21), respectively. Note that this splitting introduces some numerical 
uncertainties, since the expression in Eq. (21) is quite sensitive to small deviations in 
rpo, which in turn suffers from some arbitrariness in the exact definition (due to freeze- 
out being a process with a small but finite temporal extent rather than an immediate 
effect). 
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Figure 2: Example of the evolution of a distribution function of sterile neutrinos. One can 
clearly distinguish two momentum scales (global maximum at X\ ~ 1.5 and second, local 
maximum at £2 ~ 26). The mean comoving momentum ( x ) is located at ( x ) ~ 16.6. The 

standard deviation \J ( x 2 ) — ( x ) 2 is approximately 26, which illustrates that the mean 
value ( x ) contains practically no information. 


Estimating the two corresponding free-streaming horizons, cf. Eq. (24), 

Asir™.,) ~ (0.05,0.01) Mp C , 
bs”*“,«:r“) ~ (0.7,0.1) M P c, 


the left peak tends to yield cold/warm DM while the right one would indicate hot 
DM in both cases. The full distribution function corresponds to (Apg™®”ai aI , Apsjtotai) ~ 
(0.27, 0.05) Mpc, which is perfectly in between the two individual estimates and consis¬ 
tently indicates a case at the borderline of warm/hot DM. Even though the splitting into 
two distinct parts introduces some extra numerical uncertainty with respect to the values 
for the complete distribution function, these values clearly illustrate the issue with using 
the free-streaming horizon as an estimator. 

This is precisely one of the cases where more detailed studies about structure formation 
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Figure 3: Lines of correct abundance (dark colours) and of sizeable but insufficient abun¬ 
dance (faint colours) in the Cnp-Cr plane for different values of the sterile neutrino mass. 
In addition, the Tremaine-Gunn and overclosure bounds are displayed, see text for details. 


have to be performed in order to obtain a final answer. With the simple estimate of 
the free-streaming horizon alone, such a scenario should not be prematurely discarded. 
In any case, it is worthwhile noting that a one-component DM model can have two 
similarly important momentum scales in its distribution function, which may open up new 
possibilities to tackle the small-scale problems from structure formation simulations. 


6 Results and bounds 


In this section we will present our detailed results and we will also address possible bounds 
on the scenario as well as the validity of our considerations. 
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6.1 The Dark Matter abundance 


Let us first discuss the DM abundance, which is displayed in Fig. [3] (similarly to the one 
shown in Sec. [2|. In this plot, the dark coloured bands mark the regions in the Cnp-Cp 
plane where a sterile neutrino of a given mass yields an abundance in accordance with 
the 3 a range from the Planck 2015 data [21. For example, the dark red lines correspond 
to a sterile neutrino with a mass of = 1 keV. Here, the lines on the left correspond to 
the FIMP regime while the ones on the right correspond to the WIMP regime. We also 
mark, by the neighbouring fainter lines, the regions where a sizeable but not sufficiently 
large abundance is generated. 


There are two important bounds displayed in Fig. [3] Let us start with the so-called 
Tremaine-Gunn (TG) bound 67 , which is based on the idea that any collection of identical 
fermions must have a certain minimum phase space density. Applying this bound to the 
observed dwarf satellite galaxies leads to a lower bound of roughly mjy > 0.5 keV on the 
sterile neutrino mass 68 . Thus, smaller masses are ultimately forbidden by the Pauli 


exclusion principle. In our plots this bound is marked by the gray dashed line around the 
white area which in particular cuts into the parameter space for values of Cup close to 
one. It basically indicates where the relic density for tun = 0.5 keV equals the upper 3cr 
bound \2\. The second bound comes from the fact that the DM should not “overdose” the 
Universe, i.e., its energy density fraction Ddm should be smaller than one. Since in the 
figures we marked the lines of correct abundance for different masses m n, the resulting 
forbidden regions do in fact also depend on the mass of the DM particle. However, given 
that there is a model-independent lower value for the mass from the TG bound, at least 
for this smallest mass of m n = 0.5 keV the overclosure region marks an absolute bound. 
Not too surprisingly, this forbidden region is smaller than that excluded by the TG bound, 
and in our plots it is marked by the gray patches in the upper right corners. 


As already indicated in Secs. [2] and [4j depending on the exact values of the parame¬ 
ters different regimes are possible. Let us discuss a few numerical examples. Starting 
with the FIMP case, in Fig. [f] we illustrate an example point (CnpTr) = (10 _1 ,10 _1 ) 
which corresponds to this regime. Taking again the benchmark value ms = 1 TeV and 
using g* = 106.5 at high temperatures, the effective couplings translate into ( y , A) = 
(8.4 x 10~ 8 , 2.6 x HT 7 ) on the Lagrangian level. On the left panel, we depict the evolu¬ 
tion of the sterile neutrino distribution function /^(r) with the time parameter r. As one 
can see, most of the abundance is produced around the time r ~ 1. Soon afterwards the 
production ceases such that even for very late times, the distribution hardly changes (as 
soon as r ~ 10, the distribution is practically identical to the final one). The distribution 
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Figure 4: The evolutions of the distribution function (left) and of the abundances (right) 
for a point corresponding to the scalar being a FIMP. 


exhibits a clear peak whose maximal value is very close to the mean momentum over 
temperature, (x) ~ (p/T) ~ 2.67. This means that this sterile neutrino distribution is 


colder than a thermally produced one, for which this number would be equal to 3.15 35 


However, this point nevertheless turns out to be in the hot DM region, cf. Sec. 6.2 This 


is mainly due to the fact that this point in parameter space requires a mass of the sterile 
neutrino of ps 2keV in order to fulfil the relic abundance constraint. This low mass 
in turn leads to a long time of highly relativistic free-streaming. 


The evolution of the abundances hi, i.e., the integrals over the distribution function 
divided by T 3 , is depicted on the right panel of Fig. |4j We display both the abundances 
hs of the scalar and hjsr of the sterile neutrino. Starting with the scalar (solid black line) 
we can see that, as to be expected from a generic FIMP, the abundance of the scalar is 
gradually built up with increasing time parameter r. However, it never reaches a thermal 
abundance, as can be seen by comparing the black curve to the hypothetical one for a 
scalar in thermal equilibrium (dotted gray line). If the scalar was stable, its abundance 
would not anymore change after the freeze-in is completed, cf. dashed gray line, and it 
would in practice act as some type of DM. However, given that the scalar is unstable, 
once it does decay its abundance decreases and instead a sizeable abundance hjv of sterile 
neutrinos is built up (red solid line). Indeed, because of each scalar decaying into two 
sterile neutrinos, the final abundance of sterile neutrinos is exactly twice the one of the 
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would-be-stable scalar for late times (numerically we obtain njsr{r — 250)/ns(Cr = 0, r = 
250) ~ 2.03, in excellent agreement with the expectation). Furthermore, we can use 
Eq. (18) to cross-check our numerics, and both results agree nearly perfectly with each 

in this case. Note that this value is not a measure 


other, within a deviation of only 2. 
of the quality of our numerical methods since we do not know a priori in which part of the 
parameter space the analytical results approximate the exact result to a desired accuracy. 
We expect the deviation to become smaller as Cr further decreases. In fact, on the edge of 
our parameter space where Cr = 10~ 4 , the analytical result is reproduced with deviations 
well below 1%. 




Figure 5: The same as Fig. |4j but for the WIMP case with decay in equilibrium. 


Let us now turn to the WIMP cases. For a large enough decay width, an equilibrated scalar 
can decay while being in thermal equilibrium. This regime is in fact only realised in a rela¬ 
tively small corner of the parameter space, but one point which is a good example for this 
behaviour is (C H p,Cr) = (10 4 ,10 -L5 ) - corresponding to (y, A) = (4.72 x 10~ 8 , 8.3 x 1CT 5 ) 
for our standard reference values of ms = 1 TeV and g * = 106.5 - as depicted in Fig. [5j 
Glancing at the distribution function (left panel) first, the evolution appears to be rel¬ 
atively similar to the FIMP example just discussed - although the distribution looks 
slightly flatter at early times. This distribution also seems to be slightly colder than a 
thermal one, as to be expected [3lJ,|38,3^j, 47,48]. However, also this point will turn out 
to correspond to hot DM. 

The evolution of the abundances reveals the difference to the FIMP case, cf. right panel of 


24 


















Fig. [5] Here, it is clearly visible that the scalar (solid black line), although starting with a 
vanishing initial abundance, equilibrates very quickly and then follows the thermal curve 
(dotted gray line) |^] During this time, the scalar is highly relativistic and thus its decay 
is in fact not very efficient. However, given that its abundance is thermal and thus very 
large, the few occasional decays are sufficient to gradually built up a sizeable abundance 
of sterile neutrinos]^] The scalar remains in equilibrium for a relatively long time, until 
r ~ 5, and if it was stable it would just resemble the generic behaviour for frozen-out 
WIMPs (cf. dotted gray curve). However, given that the scalar decays relatively quickly, 
the frozen-out abundance does not survive and is converted into sterile neutrinos (solid 
red line). 


We can again compare the abundances for late times, which in this case gives hjv(r = 
250)/hg(Cr = 0,r = 250) ~ 327.6. This is vastly different from the previous result, but 
this behaviour is easy to understand: as long as the scalar is in equilibrium, it may decay 
more or less arbitrarily fast without its abundance being affected, because it is constantly 
re-generated by the thermal plasma. This happens very efficiently, so that a very large 
number of sterile neutrinos is produced while the scalar is still in thermal equilibrium. Of 
course, for the frozen-out abundance alone, the factor of two would again be present - but 
this part makes up only about 0.6% of the final sterile neutrino abundance. Hence, this 
case indeed corresponds very well to the limit of only having scalar decays in equilibrium. 
Using the approximation obtained in Eq. (20), we indeed obtain a final abundance which 
agrees with the numerical result within 4.5%. 


When we go down to smaller values of C r, we will reach the intermediate WIMP regime 
where some scalars decay in equilibrium while another non-negligible fraction decays only 
after freeze-out, thereby contributing to the final sterile neutrino abundance in simi¬ 
lar amounts. The example displayed in Fig. [6] is (C H p,Cr) = (10 4 ,10 -4 ), which corre¬ 
sponds to the bottom right corner of the parameter space considered, and to ( y , A) = 
(2.7 x 10~ 8 9 , 8.3 x 1CT 5 ) for the reference values for mg and g *. Having a look at the 
distribution function first, see left panel, one can see that a twin peak structure is visible 
-just as seen in the example for the free-streaming horizon failing, cf. Fig. [2] Looking at 
the evolution with the time parameter r, it is visible that for early times the left (lower 
momentum) peak gradually builds up, while the right (higher momentum) peaks is only 
generated much later. This behaviour already suggests that the left peak arises from 
decays in equilibrium while the right peak is generated by late decays with the scalar S 
already having frozen out and thus being out of equilibrium. This notion is supported 

8 We have explicitly checked that this happens independently of the initial abundance, as it should. 

9 In fact, one could equally well interpret this case simply as the sterile neutrino itself being a FIMP. 
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Figure 6: The same as Fig. |tJ, but for the intermediate WIMP case. 


by the late peak being the one corresponding to higher momenta: the freeze-out of the 
scalar happens at a temperature close to its mass, while the decay is a gradual process 
which takes some time to happen after the scalar has become non-relativistic. Thus, some 
energy is “stored” in the scalar mass and, once the scalars decay, the characteristic energy 
scale of the resulting sterile neutrinos will be of the order of the scalar mass - which may 
be considerably larger than the temperature of the Universe at that time. 


Glancing at the right panel we can see a matching behaviour in the abundances. While the 
scalar is in equilibrium, it gradually builds up a sizeable abundance of sterile neutrinos. 
However, this process ceases to be efficient once the scalar turns non-relativistic, thereby 
dropping in abundance and ultimately freezing out. Although the scalar abundance is 
much smaller than it was during the equilibrium time, now that the particles are non- 
relativistic the decays become efficient and completely translate the abundance of frozen- 
out scalars into sterile neutrinos, where the particle number is again doubled (but only 
for the frozen-out part). 


Finally, we have in Sec. |4.2| also discussed the case where the scalar practically fully decays 
out of equilibrium. On the other hand, given that the case just discussed was already 
located at the very edge of the parameter space, this final case does not seem to make 
sense at all. However, this does not mean that the decay solely out of equilibrium does not 
exist, but only that it is not very accurate to treat it under the assumption g* ~ const., 
as we will demonstrate in Sec. 6.5 We still want to briefly discuss a toy example of this 
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Figure 7: The same as Fig. [TJ but for the Out-of-equilibrium WIMP case. 

case, for parameter values of Chp = 10 4 and Cr = 10 -5 , which is depicted in Fig. [7] and 
which corresponds to Lagrangian level couplings of (y, A) = (2.4 x 10 10 , 8.3 x 1CU 5 ) in 
our benchmark scenario. 

Even though the double logarithmic scale in both panels might be misleading, the main 
contribution now comes from the late decay of frozen-out scalars. This can be seen easily 
when considering the sterile neutrino abundance hjv in the right panel. While the scalar is 
in equilibrium, a small abundance of sterile neutrinos builds up until a plateau is reached 
for r ~ 10. Subsequently, the decay sets in and the relic abundance of scalars is converted 
into sterile neutrinos at high momenta. The final abundance exceeds the value of the 
intermediate plateau by more than a factor of 20. Also, the average momentum ( x) is 
strongly dominated by the high momenta, just as expected. Since the production during 
equilibrium is proportional to Chp, we would have to lower this parameter by several 
orders of magnitude to make the Erst peak vanish in a log-log plot. As already argued, 
such a case would be far beyond the validity of our assumption of a constant g * during 
production and will hence not be considered in this work. 

Up to now, we have mainly been concerned with the DM abundance, but we have not yet 
shown whether the DM produced is in accordance with structure formation. A discussion 
of this point will be given in the next subsection. 
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6.2 Results for the free-streaming horizon 


As discussed in Sec. 5.1 we present the free-streaming scale 1) as calculated by our nu¬ 
merical approach, fully taking into account the evolution of a(t) (cf. Appendix [Bj, and 2) 
following the estimates put forward in 39, Eq. (20)]. In Fig. [8j we display again the bands 
in the Cup-Cr plane reproducing the correct relic abundance (just as in Figs. [I] an d§ , but 
this time colour-coding whether the sterile neutrino is hot (red), warm (purple), or cold 


(blue) according to the definitions in Sec. 5.1 The scalar mass is chosen to be ms = 1 TeV, 
however, the results depend only mildly on the mass of the scalar 39,40,47 . This is true 


since in our computation the only dependence of the scalar mass enters through the effec¬ 
tive number of d.o.f. which are a function of physical temperature. The strongest physical 
dependence on the mass of the scalar is still implicit in the definition of the parameters 
Chp and Cp. If we lowered the scalar mass to some hundred GeV, the result would be 
altered only by a few percent. 

In Fig. [8] it is clearly visible that our numerical results disfavour more of the parameter 
space than the analytical estimates do. In particular, everything but the FIMP case is un¬ 
der tension in that analysis. However, we want to emphasise once more that these results 
should be interpreted with care. First of all, the numerical approach also suffers from 
(mainly systematical) uncertainties: the simplified truncation of the time-temperature 
relation into two distinct regimes (either purely radiation dominated or purely matter 
dominated) will differ from the exact time-temperature relation, see Appendix [B] for de¬ 
tails. Second, as we will see in Sec. 6.5, in the region where the sterile neutrinos are 
particularly hot (i.e., for small decay constants), the approximation of a constant number 
of d.o.f. during production becomes less reliable, which also affects the calculation of Afs- 
Third, as discussed in Sec. 5.1, the free-streaming horizon - only taking into account 
average properties of the spectrum - cannot capture all features of structure formation 
and can hence only serve as an indication. More detailed analyses can be done using the 
so-called transfer function, i.e. the square root of the linear power spectrum of matter per¬ 
turbations, which can in turn be constrained by data from the Lyman-a measurements. 


Such studies have been performed in 31 for sterile neutrino dark matter produced from 


the DW-mechanism, from the SF-mechanism, and from a simplified version of scalar 
decay, using the Boltzmann solver CLASS 69 to obtain the transfer functions from an 


extended Press-Schechter approach. A similar study taking into account the subtleties of 
sterile neutrino DM from scalar decay as discussed in this paper is subject of on-going 
work 50 . First indications of this analysis seem to confirm the conclusions drawn from 


Fig. 8a rather than those of Fig. 8b, which puts the scenario of freezing-out scalars under 


severe tension and might also indicate a comparatively large lower mass limit of the sterile 
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(a) Numerical results for the free-streaming horizon. 



(b) Same as Fig. 8a but using the analytical estimate. 


Figure 8: Comparison of numerical results and analytical estimates of the free-streaming 
horizon, i.e., with and without taking into account the evolution of a(t) in the computation 
of the integral for Aps- Note that the strong dependence of these results on ms is hidden 
in the definition of Chp and Cr- However, there is a residual weak dependence on ms since 
the time-temperature relation used to calculate the free-streaming horizon is sensitive to 
the absolute temperature scale (as opposed to Chp and Cr). Hence we explicitly state 
the benchmark value of ms = ITeV even though the result will not change dramatically 
for scalar masses varying by a factor of a few. For this reason and for the sake of being 
able to compare to the other plots more easily, we also renounce labelling the axes with 
Lagrangian level couplings. 






























































neutrino of roughly 20 keV for the case of freeze-in, which might be an interesting finding 
in particular in the context of the tentative 3.5-keV line 19,20 . However, in order to 


give a definitive answer, we will extend the current study 49,50 , as already discussed in 
the introduction. 


6.3 The dark radiation bound 


Somewhat related to the bound from structure formation is the bound on the effective 
number of neutrinos, N e g. In the SM, this number is equal to is 3.046, the small deviation 
from 3 arising due to the effects of the reheating at e + e~ decoupling 70 . However, if the 


sterile neutrinos in our setting are too hot, they effectively act as radiation and could in 
that case also contribute to the deviation AN e g of N e s from its standard value. 

We can calculate the contribution of the sterile neutrinos to AiV e g by comparing the 
kinetic part of their energy density to the energy density °f a perfectly relativistic 
(i.e. massless) fermionic species at the same temperature in equilibrium 


AlV efr (T) = 


p — nm N 

^erm 
therm 


2 P 


60 

77T 4 


m N 

T 


dxx 2 \Jl+(^-) -1 \f N (x,T). 


m N /T 


(26) 

The factor of 2 in the denominator of the first term is due to the fact that our distribution 
function already contains both particle and antiparticle while N e g is constructed in a way 
to reproduce the number of families, i.e. 3 (up to the aforementioned small corrections) 
and not a value of 6 . Note also, that the factor ( T/T v ) 4 accounts for the fact that, once 
the reaction e + e _ -H- 77 freezes out, the photons get reheated while the neutrinos have 
already decoupled from the plasma, such that no energy is transferred to neutrinos by the 
annihilation of the electron-positron pairs. This is also true for sterile neutrinos, however, 
the temperature of their distribution (if one can define this quantity at all, given that the 
distribution may be highly non-thermal) is implicitly contained in the quantity /at(x,T). 
Yet, the relative reheating of the photons compared to sterile neutrinos is the same as that 
compared to active neutrinos, at least if the small corrections due to weak interactions are 


neglected. Thus we can simply use the factor (T/T u ) in Eq. (26), where T v is indeed the 


10 Note that this is slightly different to the standard case found in the literature, where the dark 


radiation component is typically highly relativistic, see e.g. 71 72 , whereas in our case we do not know a 


priori whether this is the case and have to take this subtlety into account by subtracting the rest energy 
which is negligibly small in the highly relativistic limit. Alternatively, one can estimate the contribution 


of non-thermal DM to dark radiation by using the Lorentz factor 73 . 
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temperature of the active neutrinos. In this case, the value of ( T/T v ) 4 rises from unity 
to (11/4) 4 ^ 3 3.85. Hence we include the factor of 3.85 for late times after electron- 

positron-annihilation while we drop it for temperatures above about 1 MeV. 


Using Eq. (26), we can - for every point (Cnp,Cr) - simply fix the sterile neutrino mass 


such as to reproduce the correct relic abundance and then calculate AN e s at any given 
temperature T. Again, all the information lies in the distribution function fw(x, T), where 
x = p/T: the higher the distribution peaks at large momenta, the bigger the contribution 
to the extra radiation will be. If, on the other hand, the sterile neutrino abundance was 


tiny, this would also be reflected in f^(x,T) and the integral in Eq. (26) would yield a 
vanishingly small result. 

In general, we have information on A N e g from two different epochs in the history of the 
Universe: during big bang nucleosynthesis (BBN), at 2~bbn = 4 Me Vp] the formation rate 
of nuclei depends on the overall expansion rate of the Universe which, in turn, depends 
on its overall radiation content 57, 77, 78 . Thus, if we do not want to spoil BBN, we 


have to respect an upper bound on the amount of extra radiation at BBN. While the 
Particle Data Group still cites a relatively ole 
newer versions exist: AiV e g BN < 1@95% C.L. 


Particle Data Group still cites a relatively old bound, AiV e BBN < 1.5@95% C.L. 179,80 

AiV BBN < 0.93@95% C.L 


81 


and 

A1V BBN < 0.85@95% C.L. 83 . We have in our plots adopted the most stringent limit, to 


illustrate that not even the strongest constraint does influence our results in a significant 
way. A seemingly more stringent constraint can be obtained from the measurement of 
the cosmic microwave background (CMB), which decouples at Tqmb ~ 0.26eV 84 


since 


the CMB spectrum also depends on the expansion rate of the Universe and thus on the 
radiation content 85 . The bound for that time is A:V/ MB < 0.32@95% C.L. |2|. 


In our analysis, we take into account both bounds and calculate AA eff (T B bn) as well as 
A N e s (Tcmb)- Although the BBN bound on A N e g appears to be weaker than the one from 
the CMB measurement, one has to take into account that BBN happens much earlier than 
the CMB decoupling. Thus, given that the sterile neutrinos produced from scalar decays 
cool down as time goes by, it may very well be that their contribution to the radiation 
content at BBN is much more significant than later on (and, indeed, this will turn out to 
be the case here). Such settings with a type of dark radiation that contributes differently 
at BBN time than later on are known in other contexts, too (see, e.g., 


Refs. 


)■ 


n The beginning of BBN happens at a temperature of a few MeV 74 , and we take 4 MeV as example 

Given that the sterile 


which is known to “reset” conditions to how they were prior to BBN 75 76 


neutrinos keep cooling down until the end of BBN and even further, taking such an early temperature 
corresponds to some extend to a pretty conservative limit. However, we also have to take into account 
that the main temperature dependence is factored out of AN e g per definition, so that Tbbn = 4 MeV is 
in fact not much more conservative than taking a value of 1 MeV, or similar. 
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Deviation of the effective number of neutrinos from the SM 
Le at BBN and resulting excluded region (red patches). 



(b) Same as Fig. 9a at the temperature of CMB decoupling. 


Figure 9: Deviation of the effective number of neutrinos from its SM value of 3.046. 



















Some example contour lines for AN e g are displayed in the Cnp-Cr plane in Fig. [9j both 
for BBN (upper panel) and CMB (lower panel). As a guide for the eye, we have again 
displayed the lines of correct abundance, cf. Fig. [3j since we fix the mass of the sterile 
neutrino entering the computation of A N e g via Eq. (26) by the constraint of reproducing 
the observed relic abundance. As can be seen, at BBN there could be a non-negligible 
contribution of the sterile neutrinos to AA e ff (Xbbn), and there is even a small excluded 
region which would violate the bound (marked by the red areas at the bottom centre of 
the plot). This region of a too large contribution arises only for very small decay widths 
Cp, i.e., the scalars must be very long-lived and inject highly energetic sterile neutrinos 
into the Universe at relatively late times. However, we would in any case exclude this 
region from any serious consideration because it would, trivially, fall far into the region 
where the DM is hot in any case, cf. Fig. [8] At CMB, on the other hand, there is not 
even a serious constraint left since the sterile neutrinos have cooled down by then and 
only have comparatively small momenta. 


Thus, even though there is in principle a contribution of the sterile neutrinos to AN e g, no 
strong constraint arises from it and there is no threat to our production mechanism. 


6.4 Other bounds and constraints 


In this section, we will discuss the two remaining conditions which may affect the DM 
production mechanism proposed. It turns out that both of them are no actual problems: 
the first one would only affect regions so far away from the interesting part of the pa¬ 
rameter space that they do not play a role in practice. The second problem is more of 
a “theoretical” nature, i.e., while it may be important to take into account, it can be 
easily circumvented in concrete settings and may only appear to be problematic if the 
Lagrangian presented in Eq. (|TJ) was viewed as “theory of everything” valid up to the 
Planck scale. However, for completeness we would like to at least briefly mention these 
points. 


In general, collider bounds could also restrict the parameter space of our model, since 
after all the Higgs portal coupling A could be used to produce two singlet scalars from 
two SM-like Higgses. Using the limits on the mixing between a scalar singlet and the 
Higgs boson as proposed in 90 , the bounds on A are far above even the largest value of 
Chp we show in our plots. This holds true even if the VEV of the scalar singlet is larger 
than its mass by two orders of magnitude. This behaviour can be intuitively understood 
by taking into consideration that even a value of Chp = 10 4 corresponds to rather small 
couplings A, due to the large factor M 0 /ms involved in its definition. Thus, in practice, 
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we do not need to be bothered by any current bounds from colliders. However, at least 
in principle, there is an upper bound on Chp which may become important if our study 
was extended to considerably larger values of the Higgs portal coupling. This conclusion 


still holds when confronted with updated analyses 91 


There is an orthogonal problem which is related to the symmetry breaking resulting from 
the scalar potential in Eq. ([2]). As we had already explained, the shape of this potential 
is determined by an underlying discrete Z 4 symmetry. While we only use this symmetry 
for simplicity and we could even skip it without too drastic consequences, at least in 
principle it would be broken by a VEV (S) of the scalar field. Since the potential has two 
perfectly degenerate minima, different parts of the Universe could then have their ground 
states in different minima, thereby leading to so-called domain walls 92 . The existence 


of such objects would have considerably changed the history of the Universe and they are 
hence problematic. However, there are many possible solutions to this problem, since the 
slightest difference in energies of the two vacua could be generated by all kinds of physics, 
in which case the walls would decay exponentially quickly 93-97 . This can also happen 
in the case at hand 
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6.5 Assessing the validity of approximating the relativistic d.o.f. 
as constant 

As stated at the very beginning, we followed the assumption that g*,g*s are constant 
during the DM production process. This assumption impacts on the form of the spectrum, 
cf. Eq. (J7|), and thereby also on the implications for structure formation. Even in the 
analytical estimate of the free-streaming horizon, it makes a difference if the dilution 
factor is calculated from a starting point of g* = 106.75 or at some lower value. In 
order to assess this approximation, we have for every point in the Cnp-Cp plane computed 
the temperature when the production of sterile neutrinos is completed (i.e., when the 
abundances surpasses 95% of its maximum value) and we plot the comparison to the SM 
number of d.o.f. of the primordial plasma at that temperature. To this end, we again 
assume a scalar mass of 1 TeV. 

In our case, the extension of the SM will never contribute more than l + 2x7/8 = 22/8 = 
2.75 (i.e., less than 3% of the SM value of 106.75) units to the count of d.o.f. since this 
would be the maximum possible contribution of scalars and sterile neutrinos both being 
present with a thermal abundance and relativistic velocities. Since we expect the SM 
d.o.f. to be much larger during the time of production, we can interpret the further d.o.f. 
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Figure 10: Deviation of the number of d.o.f. at the time of DM production from the 
maximum value of g * = 106.75. This deviation can be interpreted as a check of the 
goodness of the approximation of g* = const, during production, which was used to 
simplify Eq. ([7]). 


as small additional perturbation in the background of the SM. 

Hence the approximation of constant d.o.f. can be assessed by checking the evolution of 
the SM background. If the number of SM d.o.f. at end of production is still close to 


the maximum value of g* = 106.75, the approximation can be seen as adequate. Fig. 10 


shows the deviation of the number of d.o.f. at time of production from the maximum 
value of 106.75. Comparing to Figs. [8j in most of the interesting parameter space our 
approximation is not too bad. In fact, for the cooler regions in the FIMP case, the 
approximation is even excellent. Only for the relatively hot regions the estimate of the 
deviation compared to the exact treatment is larger than 10%, but this region is in any 
case not the favoured one. Since, after all, the Boltzmann treatment of DM production 
in any case cannot be expected to yield sub-percent accuracy 98 , this means that the 
approximation is in fact not too bad. 


35 















7 Conclusions & outlook 


We have presented a fully comprehensive study of the production of keV sterile neutrino 
Dark Matter in the early Universe by singlet scalar decays. The current paper lays the 
foundation for several follow-up considerations. Aiming at a clear overview first, we have 
applied some approximations that enabled us to present analytical results in addition 
to a detailed numerical computation, which we used to further back up certain limiting 
cases. Based on these initial considerations, we have derived the system of Boltzmann 
equations to be solved at the level of momentum distribution functions, and we have 
furthermore introduced a very efficient parametrisation to do so. After presenting our 
analytical results and discussing some aspects related to cosmological structure formation, 
we turned to present our numerical results. Our numerical solutions for the distribution 
functions have not only provided a comprehensive picture of where the observed Dark 
Matter abundance can be obtained in the large parameter space investigated, but they 
have also allowed us to account for all bounds applicable. This is the first time that 
such detailed results have been obtained for the production mechanism at hand, and we 
have shown how to fully exploit the information contained in the distributions. We have 
in particular found situations in which highly non-trivial distribution functions featuring 
more than one momentum scale can result from the simple decay mechanism presented, 
which could be very interesting for cosmological structure formation. In such cases, the 
simple minded estimate of structure formation properties using only the free-streaming 
horizon will fail, as we have explained by an illustrative example. While preliminary 
results of our on-going work (beyond what is presented in this paper) indicate that the 
numerical estimate of the free-streaming horizon yields a rather comprehensive picture 
despite the difficulties involved, we nevertheless have to postpone a definitive conclusion 
to the pending results. 


In general, we have not only found very good agreement between our analytical and 
numerical results, but we have also shown that the assumptions applied (in particular 
the effective number of degrees of freedom being constant during Dark Matter production 
and the singlet scalar having a mass larger than that of the Higgs boson) are good in 
a significant fraction of the parameter space, although of course certain regimes exist 
in which the error introduced gets unacceptably large. Thus, the natural next step will 
be to extend our (numerical) considerations to the regimes in which the approximations 
applied are not valid |49j, which will in particular allow us to treat considerably smaller 
singlet scalar masses. We will furthermore extend our studies to investigate the detailed 
implications for structure formation 50 , which we have clearly shown to require a more 
advanced machinery than the free-streaming horizon. Ultimately, we aim to provide a fully 
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comprehensive study of keV sterile neutrino Dark Matter production by scalar decays, so 
that this production mechanism can be put to the acid test to determine how good an 
alternative to resonant production it truly is. 
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A Appendix: Details on the kinetic equations 


The collision term for a species <f> in contact with other species via some interaction of 
the generic type $ + a + 6 + ... a + /3 + ... is given by 


C[U\ = — l — I [dP a dP b ...dP Q dPp... x (2 tt) 4 <5 (4) (p* +p a +p b + ... -p a -pp- ■ ■■) x |Af| 2 

X [faf/3-h (1 ± fa) (1 ± fb) ••• - fafb (1 ± fa) (1 ± //}) ••• (1 ± /*)]] . 

(A-l) 


Some remarks about Eq. (A-l) are in order: 


1. The quantity E Px denotes the energy of particle x and is hence given by E x = 

VpI + 


ml. 


2. The internal degrees of freedom of a species are denoted by g x . 

3. We have introduced a symbol for the invariant phase-space element: 

d 3 p x 


d-P x 9x 


2 E Px ( 2 t rf 


(A- 2 ) 


4. The plus signs apply in the case of bosons and the minus signs in the case of fermions. 
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5. The squared matrix element is defined following the convention in 84 chapter 5], 
i.e. the squared matrix element contains all relevant symmetry factors and averages 
over both the initial and the final state spins. 


A.l Kinetic equation for the scalar 

In this appendix, we derive the kinematic functions T and Q in their exact form. We also il¬ 
lustrate some pedagogical steps to ease the derivation of the relevant collision terms. 

Let us start by constructing the collision terms in the variables x and r, cf. Eq. © , under 


the further simplification that we approximate all factors in Eq. (A-l) arising from the 


final states by unity, i.e. we neglect the bosonic/fermionic nature of the particles, which 
is a good approximation provided that their energy is much larger than the temperature 
of the plasma. There are several processes that can contribute to the production of the 
singlet scalars in the early Universe. 

We start with the collision term describing the production of a pair of scalars SS from a 
pair of SM-like Higgses hh : 


Chh —>ss (?) 


2 E n 



d 3 q' d 3 p 


d 3 p' 


(2vr) 3 2 E q/ (2vrf 2 E p (2t rf 2 E p 


4A 2 (2tt) 


(A-3) 


x S (E, + Eg' -E p - E P ') 6 < 3 > (7+ <?- p-ft) C (q) f? (,') 


In Eq. (A-3), the momenta q and q' (p and p') belong to the Higgs bosons in the initial 


state (to the scalars S in the final state). In Eq. (A-3), we have explicitly inserted the 
squared matrix element \Ai\~ = 4A 2 . 


Using an argument based on detailed balance 99 , we can make the following replacement 


in Eq. (A-3): 


/r (?) /r (?o = r s q (?) r s q (p) 


(A-4) 


Note that this can be shown quite easily in an explicit way in the case of a Maxwell- 
Boltzmann approximation that we are using, exploiting only energy conservation. 

Accordingly, we can also use the explicit form of a Maxwellian distribution to simplify 


Eq. (A-3) further. Integrating out the phase space in p and p ', we obtain 


Chh ->-ss (?) — 


4A 2 1 
8n 2 E n 


d 3 </ / ( E q + E q i) 2 — qq' cos 9 — 4 m\ 


(27t) 3 2E q 


(. E q + E q 'Y - qq' cos 9 


H e -{E q +E ql )/T _ ( A _5) 
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In order to arrive at Eq. (A-5) by integrating out the phase spaces in p and p' in Eq. (A-3), 


we have used the standard phase space integral: 


d 3 P d V x4c( 4 )/ , / a 

-q-q-( 2 vr) 0 y ' (q + q “ p — p ) — 

(2vr) 3 2 E p (2vr) 3 2 E p , 


In our case, the centre-of-mass velocity is given by 


dO 


f 2^cm\ 
47T 87T V E cm J 


(A-6) 


Pcm _ 1 ( Eg + E q i) 2 - qq' cos 9 - 4mf 

-^cm 2 y (. E q + E q /) 2 — gg' cos 6 1 


(A-7) 


With this it is straightforward to arrive at Eq. (A-5) from where the definition of T can 


directly be read after changing to the variables x and r: 


T (x,r,£) 


= 2tt dxx 2 / dcos#- 


3 — Vx 2 +r 2 


-1 


(\/x 2 + r 2 + y/x 2 + r 2 ) 2 — xx cos 6 1 — 4£ 2 r 2 


x/^Tr 2 ^ (a /£ 2 + r 2 + \/x 2 + r 2 ) 2 — 


xx cos 9 


(A- 8 ) 


where £ = m-h/ms (cf. Sec. [4]). The maximum allowed value for the cosine of 9 comes 
from the trivial constraint that the argument in the square root must be non-negative. 
This leads to 


a™* = mm 


1 , max 


- 1 , 


(\/x 2 + r 2 + yjx 2 T r 2 ) 2 — 4£ 2 r 2 
xx 


(A-9) 


For scalar masses much larger than the Higgs mass, i.e. £ C l,we can simplify Eq. (|A- 8 |): 

T (x, r, £ « 1) = AnK l (r) , 

where K\ is the first modified Bessel function of second kind. 


(A-10) 


With this result at hand, the first part of the kinetic equation for the scalar, cf. Eq. 
reads: 


dr 


, - = M o 1 A 2 

dr dT hh ^ ss ms y / x 2 _|_ r 2 647T 4 ^ 


(^—Vx 2 + r 2 j J 7 (x,r) 


(A-ll) 
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For the process of a pair of scalars annihilating into a pair of Higgs bosons, the collision 
term reads: 


^SS—thh (?) — 


2 E n 


d V 


d 3 p 


dV 4A 


(2vr) 3 2£g, (2 tt) 3 2£ p (2 tt) 3 2£ p / 8vr 


X (2vr) 4 <5 (4) {q + q -p~ p) fs (q) fs W) 


fs(q ) 4A 2 f d 3 q' 

2 E q 16t tJ (2n) 3 2E q , 


(A - 12) 

(F/ g + E q i) -qq'cosO 


Again, we can directly infer the definition of Q as in Eq. (A-14): 


Q{x,r,i) 


= 27r / dx x 


OO dmax 

2 


dcos#- 


-l 


(a/x 2 + r 2 + a/x 2 + r 2 ) ^ — xx cos 9 — 4£ 2 r 2 


\/£ 2 + r 2 ^ (\/x 2 + r 2 + \/x 2 + r 2 )' 


xx cos 0 


(A-13) 


Note that the only difference between Eqs. (A-13) and (A-8) is the exponential factor in 


the integral, stemming from the equilibrium distribution in Eq. (A-3). The entity a max is 


again defined as in Eq. (A-9). 


Thus, the second part of the kinetic equation of the scalar (again cf. Eq. ^8j)) reads: 

rSS ^ hh ( r , x) = dT dt r g _ _M 0 1 A 2 

dr dT ss ^ hh ms ^j x 2 _|_ r 2 647T 4 ' 


dfl 


dr 


;fs(x,r) / d 3 xf s (x,r)g (x,r). 


(A-14) 

Note that the collision term in Eq. (A-14) contains the actual distribution function of the 


scalar on the right-hand side, yielding an integro-differential equation for the distribution 
function fg we are interested in. 

The term describing the decay of a scalar into a pair of sterile neutrinos is constructed 
as: 


Cs^-NN (q) ~ 


2 E„ 


2 d 3 p 2d 3 p' 1 9 

-~y 2 E p E p . 


(27t) 3 2E p (27t) 3 2E p / 2 l 


1 - 


P * P' 
-f—J p J—/pf 


= -fs(q) r^, 

E q 


(27r) 4 (f (4) (q - p - p') fs (q) 

(A-15) 
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2 

where we have explicitly inserted the decay width T = . This collision term can 

be interpreted intuitively: it is clear that the rate at which scalars are depleted due to 
decays must be proportional to the time-dilated decay width which involves an 

additional boost factor owing to the physical momentum q of the scalar, and that it must 
be proportional to the amount of scalars present at that particular momentum, i.e. to 

fs(q)- 

With this, the third part of the kinetic equation of the scalar, is then given by: 


dfg^ NN (r, x) 
dr 


— — C 5 - M ° 1 r 2 F A (x r) 
dr d T C ' S ^ NN ~ m s ^ m / s V’ r > 


(A-16) 


A.2 Kinetic equation for the sterile neutrino 


We also want to show the most import steps to derive at the kinetic equation of the sterile 
neutrino. The source term looks like: 


Cs-^-NN — 2 X 


1 


2 d 3 p' d 3 ps 


2 EpJJ (2vr) 3 2Ep, (2vr) 3 2E ps 
where the matrix element for the decay reads: 


(2n) a S (E rs - E p - E?) i5 (3) (ps-p-ity 2\M\ ! f s ( ps,t) , 

(A-17) 


1 2 1 o / 1 


\M\ 2 =-y 2 p-p' = -y 2 E p E p , 


1 - 


p ■ p' 


(A-18) 


According to our conventions for symmetry factors, we average over initial and final 
states. 

Integrating out the phase space in p', we can get rid of the spatial 5-distribution in 
Eq. (A-17). Doing so, we can replace p ■ p' by p ■ ps — p 2 . The scalar product p ■ ps is 


restricted by the kinematics of the process. Using —1 < cos[<(p, p^)] < 1, where <(a,b) 
is the angle between the 3-vectors a and b, this gives the constraint 


Ps > 


p- 


mt 


4 p 


— Pmin j 


(A-19) 


which implies the lower boundary in Eq. (A-20) 
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Hence, using Eq. (|9j), the kinetic equation of the sterile neutrino is finally given by 

cS^-NN _\ „2 r°° 


dfk 


dr 


Or) = 2Q r^_ 


d£ 


x 


X z 


\J X 2 + ' 


~Js(x,r), 


(A-20) 


which we use as our master equation, cf. Eq. (16). 


B Appendix: Some background on the numerical 
calculation of the free-streaming horizon 


This appendix briefly summarises some technical details implemented in our numerics 
in order to evaluate the integral occurring in the definition of the free-streaming hori¬ 
zon, cf. Eq. (25). We treat the thermal history of the Universe in a way where there 
are two distinct periods of interest, namely a purely radiation-dominated one which 
is followed by an immediate turnover into complete matter domination at (T eq , t eq ) = 
(1.48 eV, 6.04 x 10 11 s). Note that these quantities stemming from our rather crude ap¬ 
proximation agree fairly well with the values that can be found in the literature. 


During radiation dominance, we can infer the time-temperature relation from Eq. (J5]) 
using the evolution of the d.o.f. as given in |6l]. A relation between the scale factor and 
the temperature can be established by solving for T in the equation of conservation of 
comoving entropy density, 


2t r 
45 


g*sT i a i = const. = so, 


(B-l) 


normalising today’s scale factor to unity. 


evolution of the d.o.f. as presented in 61 


In our numerics, we implement the numerical 


During matter dominance, integrating the Friedmann equation gives a relation between 
the cosmological time and the scale factor which reads 


t Aq T 




or 1/2 f„3/2 


r • 


(B-2) 


This can be converted into a time-temperature relation by virtue of Eq. (B-l). 

We have checked that our numerical treatment of the time-temperature relation repro¬ 
duces other known benchmark points (like the time of CMB-decoupling) to a very rea- 
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sonable accuracy. 
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